{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Segmentation: Thresholding and Edge Detection</h1>\n",
    "\n",
    "In this notebook our goal is to estimate the location and radius of spherical markers visible in a Cone-Beam CT volume.\n",
    "\n",
    "We will use two approaches:\n",
    "1. Segment the fiducial using a thresholding approach, derive the sphere's radius from the segmentation. This approach is solely based on SimpleITK.\n",
    "2. Localize the fiducial's edges using the Canny edge detector and then fit a sphere to these edges using a least squares approach. This approach is a combination of SimpleITK and scipy/numpy.\n",
    "\n",
    "Note that all of the operations, filtering and computations, are natively in 3D. This is the \"magic\" of ITK and SimpleITK at work.\n",
    "\n",
    "The practical need for localizing spherical fiducials in CBCT images and additional algorithmic details are described in:\n",
    "Z. Yaniv, \"Localizing spherical fiducials in C-arm based cone-beam CT\", *Med. Phys.*, Vol. 36(11), pp. 4957-4966."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "\n",
    "%run update_path_to_download_script\n",
    "from downloaddata import fetch_data as fdata\n",
    "\n",
    "%matplotlib notebook\n",
    "import gui\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "from scipy import linalg\n",
    "\n",
    "from ipywidgets import interact, fixed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the volume, it contains two spheres. You can either identify the regions of interest (ROIs) yourself or use the predefined rectangular regions of interest specified below ((min_x,max_x), (min_y, max_y), (min_z, max_z)).\n",
    "\n",
    "To evaluate the sensitivity of the algorithms to the image content (varying size and shape of the ROI) you should identify the ROIs yourself. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "spherical_fiducials_image = sitk.ReadImage(fdata(\"spherical_fiducials.mha\"))\n",
    "roi_list = [((280,320), (65,90), (8, 30)),\n",
    "            ((200,240), (65,100), (15, 40))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a GUI to specify a region of interest. The GUI below allows you to specify a box shaped ROI. Draw a rectangle on the image (move and resize it) and specify the z range of the box using the range slider. You can then view the ROI overlaid onto the slices using the slice slider. The toolbar on the bottom of the figure allows you to zoom and pan. In zoom/pan mode the rectangle interaction is disabled. Once you exit zoom/pan mode (click the button again) you can specify a rectangle and interact with it.\n",
    "\n",
    "We already specify two ROIs containing the two spheres found in the data (second row below). \n",
    "\n",
    "To evaluate the sensitivity of the two approaches used in this notebook you should select the ROI on your own and see how the different sizes effect the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roi_acquisition_interface = gui.ROIDataAquisition(spherical_fiducials_image)\n",
    "roi_acquisition_interface.set_rois(roi_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the user specified ROIs and select one of them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "specified_rois = roi_acquisition_interface.get_rois()\n",
    "# select the one ROI we will work on\n",
    "ROI_INDEX = 0\n",
    "\n",
    "roi = specified_rois[ROI_INDEX]\n",
    "mask_value = 255\n",
    "\n",
    "mask = sitk.Image(spherical_fiducials_image.GetSize(), sitk.sitkUInt8)\n",
    "mask.CopyInformation(spherical_fiducials_image)\n",
    "for x in range(roi[0][0], roi[0][1]+1):\n",
    "    for y in range(roi[1][0], roi[1][1]+1):        \n",
    "        for z in range(roi[2][0], roi[2][1]+1):\n",
    "            mask[x,y,z] = mask_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholding based approach\n",
    "\n",
    "To see whether this approach is appropriate we look at the histogram of intensity values inside the ROI. We know that the spheres have higher intensity values. Ideally we would have a bimodal distribution with clear separation between the sphere and background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "intensity_values = sitk.GetArrayViewFromImage(spherical_fiducials_image)\n",
    "roi_intensity_values = intensity_values[roi[2][0]:roi[2][1],\n",
    "                                        roi[1][0]:roi[1][1],\n",
    "                                        roi[0][0]:roi[0][1]].flatten()\n",
    "plt.figure()\n",
    "plt.hist(roi_intensity_values, bins=100)\n",
    "plt.title(\"Intensity Values in ROI\")\n",
    "plt.show()                                                                                               "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you identify the region of the histogram associated with the sphere?\n",
    "\n",
    "In our case it looks like we can automatically select a threshold separating the sphere from the background. We will use Otsu's method for threshold selection to segment the sphere and estimate its radius. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set pixels that are in [min_intensity,otsu_threshold] to inside_value, values above otsu_threshold are\n",
    "# set to outside_value. The sphere's have higher intensity values than the background, so they are outside.\n",
    "\n",
    "inside_value = 0\n",
    "outside_value = 255\n",
    "number_of_histogram_bins = 100\n",
    "mask_output = True\n",
    "\n",
    "labeled_result = sitk.OtsuThreshold(spherical_fiducials_image, mask, inside_value, outside_value, \n",
    "                                   number_of_histogram_bins, mask_output, mask_value)\n",
    "\n",
    "# Estimate the sphere radius from the segmented image using the LabelShapeStatisticsImageFilter.\n",
    "label_shape_analysis = sitk.LabelShapeStatisticsImageFilter()\n",
    "label_shape_analysis.SetBackgroundValue(inside_value)\n",
    "label_shape_analysis.Execute(labeled_result)\n",
    "print(\"The sphere's location is: {0:.2f}, {1:.2f}, {2:.2f}\".format(*(label_shape_analysis.GetCentroid(outside_value))))\n",
    "print(\"The sphere's radius is: {0:.2f}mm\".format(label_shape_analysis.GetEquivalentSphericalRadius(outside_value)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception thrown in SimpleITK Show:"
   },
   "outputs": [],
   "source": [
    "# Visually evaluate the results of segmentation, just to make sure. Use the zoom tool, second from the right, to \n",
    "# inspect the segmentation.\n",
    "gui.MultiImageDisplay(image_list = [sitk.LabelOverlay(sitk.Cast(sitk.IntensityWindowing(spherical_fiducials_image, windowMinimum=-32767, windowMaximum=-29611),\n",
    "                                      sitk.sitkUInt8), labeled_result, opacity=0.5)],                   \n",
    "                      title_list = ['thresholding result'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Based on your visual inspection, did the automatic threshold correctly segment the sphere or did it over/under segment it?\n",
    "\n",
    "If automatic thresholding did not provide the desired result, you can correct it by allowing the user to modify the threshold under visual inspection. Implement this approach below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Edge detection based approach\n",
    "\n",
    "In this approach we will localize the sphere's edges in 3D using SimpleITK. We then compute a least squares sphere that optimally fits the 3D points using scipy/numpy. The mathematical formulation we use is as follows:\n",
    "\n",
    "Given $m$ points in $\\mathbb{R}^n$, $m>n+1$, we want to fit them to a sphere such that\n",
    "the sum of the squared algebraic distances is minimized. The algebraic distance is:\n",
    "$$\n",
    "\\delta_i = \\mathbf{p_i}^T\\mathbf{p_i} - 2\\mathbf{p_i}^T\\mathbf{c} + \\mathbf{c}^T\\mathbf{c}-r^2\n",
    "$$\n",
    "\n",
    "The optimal sphere parameters are computed as:\n",
    "$$\n",
    "[\\mathbf{c^*},r^*] = argmin_{\\mathbf{c},r} \\Sigma _{i=1}^m \\delta _i ^2\n",
    "$$\n",
    "\n",
    "setting $k=\\mathbf{c}^T\\mathbf{c}-r^2$ we obtain the following linear equation system ($Ax=b$):\n",
    "$$\n",
    "\\left[\\begin{array}{cc}\n",
    "      -2\\mathbf{p_1}^T & 1\\\\\n",
    "      \\vdots        & \\vdots \\\\\n",
    "     -2\\mathbf{p_m}^T & 1\n",
    "     \\end{array}\n",
    "\\right]\n",
    "\\left[\\begin{array}{c}\n",
    "      \\mathbf{c}\\\\ k\n",
    "      \\end{array}\n",
    "\\right] =\n",
    "\\left[\\begin{array}{c}\n",
    "      -\\mathbf{p_1}^T\\mathbf{p_1}\\\\\n",
    "      \\vdots\\\\\n",
    "      -\\mathbf{p_m}^T\\mathbf{p_m}\n",
    "      \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "The solution of this equation system minimizes $\\Sigma _{i=1}^m \\delta _i ^2 = \\|Ax-b\\|^2$.\n",
    "\n",
    "Note that the equation system admits solutions where $k \\geq\n",
    "\\mathbf{c}^T\\mathbf{c}$. That is, we have a solution that does not\n",
    "represent a valid sphere, as $r^2<=0$. This situation can arise in\n",
    "the presence of outliers.\n",
    "\n",
    "Note that this is not the geometric distance which is what we really want to minimize and that we are assuming that there are no outliers. Both issues were addressed in the original work (\"Localizing spherical fiducials in C-arm based cone-beam CT\").\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create a cropped version of the original image.\n",
    "sub_image = spherical_fiducials_image[roi[0][0]:roi[0][1],\n",
    "                                      roi[1][0]:roi[1][1],\n",
    "                                      roi[2][0]:roi[2][1]]\n",
    "\n",
    "# Edge detection on the sub_image with appropriate thresholds and smoothing.\n",
    "edges = sitk.CannyEdgeDetection(sitk.Cast(sub_image, sitk.sitkFloat32), lowerThreshold=0.0, \n",
    "                                upperThreshold=200.0, variance = (5.0,5.0,5.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Get the 3D location of the edge points and fit a sphere to them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "edge_indexes = np.where(sitk.GetArrayViewFromImage(edges) == 1.0)\n",
    "\n",
    "# Note the reversed order of access between SimpleITK and numpy (z,y,x)\n",
    "physical_points = [edges.TransformIndexToPhysicalPoint([int(x), int(y), int(z)]) \\\n",
    "                   for z,y,x in zip(edge_indexes[0], edge_indexes[1], edge_indexes[2])]\n",
    "\n",
    "# Setup and solve linear equation system.\n",
    "A = np.ones((len(physical_points),4))\n",
    "b = np.zeros(len(physical_points))\n",
    "\n",
    "for row, point in enumerate(physical_points):\n",
    "    A[row,0:3] = -2*np.array(point)\n",
    "    b[row] = -linalg.norm(point)**2\n",
    "\n",
    "res,_,_,_ = linalg.lstsq(A,b)\n",
    "\n",
    "print(\"The sphere's location is: {0:.2f}, {1:.2f}, {2:.2f}\".format(*res[0:3]))\n",
    "print(\"The sphere's radius is: {0:.2f}mm\".format(np.sqrt(linalg.norm(res[0:3])**2 - res[3])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception thrown in SimpleITK Show:"
   },
   "outputs": [],
   "source": [
    "# Visually evaluate the results of edge detection, just to make sure. Note that because SimpleITK is working in the\n",
    "# physical world (not pixels, but mm) we can easily transfer the edges localized in the cropped image to the original.\n",
    "# Use the zoom tool, second from the right, for close inspection of the edge locations.\n",
    "\n",
    "edge_label = sitk.Image(spherical_fiducials_image.GetSize(), sitk.sitkUInt16)\n",
    "edge_label.CopyInformation(spherical_fiducials_image)\n",
    "e_label = 255\n",
    "for point in physical_points:\n",
    "    edge_label[edge_label.TransformPhysicalPointToIndex(point)] = e_label\n",
    "\n",
    "gui.MultiImageDisplay(image_list = [sitk.LabelOverlay(sitk.Cast(sitk.IntensityWindowing(spherical_fiducials_image, windowMinimum=-32767, windowMaximum=-29611),\n",
    "                                                                sitk.sitkUInt8), edge_label, opacity=0.5)],                   \n",
    "                      title_list = ['edge detection result'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## You've made it to the end of the notebook, you deserve to know the correct answer\n",
    "\n",
    "The sphere's radius is 3mm. With regard to sphere location, we don't have a the ground truth for that, so your estimate is as good as ours. "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
